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Abstract. We present here an extension of the Wang-Landau Monte Carlo 
method which allows us to get very accurate estimates of the full probability 
distributions of several observables after a quantum quench for large systems, 
whenever the relevant matrix elements are calculable, but the full exponential 
complexity of the Hilbert space would make an exhaustive enumeration impossible 
beyond very limited sizes. We apply this method to quenches of free-fermion 
models with disorder, further corroborating the fact that a Generalized Gibbs 
Ensemble fails to capture the long-time average of many-body operators when 
disorder is present. 
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1. Introduction 

A sudden quench of the Hamiltonian parameters is perhaps the simplest form of out- 
of-equilibrium dynamics that a closed quantum system can experience. Experiments 
with “virtually isolated” cold atomic species in optical lattices mm have transformed 
this seemingly theoretical dream into a rich and lively stage. Several fundamental 
issues of theoretical quantum statistical physics, like the onset of thermalization which 
is generally expected to occur for a closed quantum system after a sudden quench 
HEME or the “breakdown of thermalization” mm expected when the system is 
integrable or nearly integrable, are now of experimental relevance [8j. We refer the 
reader to a recent review 0] for an extensive introduction to such non-equilibrium 
quantum dynamics issues. 

The issue we want to tackle in this paper is the following. Suppose you perform 
a quantum quench of the Hamiltonian parameters, abruptly changing, at t = 0, from 
Hq —» H. If l^o) denotes the initial quantum state at t = 0, and |ct) the eigenstates 
of H with energy E a , the ensuing quantum dynamics would lead to averages for any 
given operator A given by: 

A(t) = (Vole^Ae-^o) = £ |c«| 2 A act + £ e« E «'- E ^c* a ,A a , a c a , (1) 

o: 

where c a = (aj'Lo) and A a / a = (c/|A|a). The first (time-independent) term in the 
previous expression dominates the long-time average of A(t), and is usually known as 
diagonal average [5) 

(i) D = ^|c a | 2 A aQ . (2) 

Oi 

To calculate it, in principle, we should take the sum over all the (many-body) 
eigenstates |a) of H — an exponentially large number of states —, calculating for 
each of them the overlap c a and the associated diagonal matrix element A aa . Luckily, 
(A) d can be calculated for many problems, notably those that can be reduced to 
quadratic fermionic problems, by circumventing in one way or another exponentially 
large sums: for instance, through a detour to time-dependent single-particle Green’s 
function and the use of Wick’s theorem, see e.g. mm- But suppose that you want 
to know more than just the diagonal average (A) D , and pretend to have information 
on the whole distribution of the values of A aa accessed after the quench mm, i-e., 

Po (A) = J2\c a \ 2 S(A-A aa ) , (3) 
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of which (A) D is just the average: (A) D = fdAp D (A) A. Here there is, evidently, a 
problem: knowing the distribution of A requires exploring the full many-body Hilbert 
space, summing over the eigenstates |a), and this exhaustive enumeration would 
restrict our calculations to exceedingly small sample sizes, although all information 
on c a and A aa might in principle be easy to calculate, or in any case accessible, for 
instance by just solving a one-body problem (hence, for much larger sizes). A similar 
problem occurs in considering, for instance, the corresponding [f] generalized Gibbs 
ensemble (GGE) average [HI 1T51 [T5j HT j fig] 

v—.. e~ ^-G GG 

(^4)gge = ^2 y A aa , (4) 

Z GGE 


where are Lagrange multipliers which constrain the mean value of each of the 


constants of motion to their t = 0 value, (\E' 0 |/ M |'I'o) = Tr 


Pgge-I/j 


> = («l4l a )> 


■Z’gge is the GGE partition function, and p GGE = e~ /Z GGE . Once again, for 

“quadratic problems” this average is rather simply calculated in terms of single-particle 
quantities, but the corresponding distribution 


A I a 

Pggv(A) = ^2 - S(A-A aa ), (5) 

a Z GGE 

requires a difficult sum over the Hilbert space. [|] 

Concerning the issue of thermalization after a quantum quench, we might indeed 
expect that, if the system is well described by a GGE ensemble, not only the mean 
values of p E (A) and p GGB (A) are equal, i.e., (A) D = f d Ap E (A)A = (A) GGE = 
f d A p GGE (A)A, but also the two distributions should be closely related; at least this 
is what a good statistical ensemble should do. 

Quite generally, we might formulate the problem as follows: how can we obtain 
information on weighted distributions (or density of states) 


Pw(A) = ^2 w a S(A - A a ) , 

a 


( 6 ) 


with positive weights w a , when both w a and A a are “easily calculated”, but the 
sum over a runs over an exponentially large “configuration space”? As discussed 
before, examples of this are the diagonal distribution, where w a = \c a \ 2 : the GGE 
distribution, where w a = e~-?G GG /Z GGE , but also the microcanonical distribution, 
where w a is a window characteristic function for the microcanonical shells, etc. A 
Monte Carlo algorithm to perform such exponentially large sums in configuration space 
seems unavoidable. We stress that this is so even if one is considering quenches in 
free-fermion models, where the relevant many-particle states |a) and matrix elements 
are easy to write down and calculate. The alternative of using exact diagonalization 
methods would put a strong limit on the size of the problem which can be studied. 

In this paper we introduce a Monte Carlo method — obtained by a rather natural 
extension of the Wang-Landau algorithm (WLA) [HI [20 [Hi — which will allow us to 
compute weighted distributions of the form of Eq. ([h]). The Wang-Landau algorithm, 


| The generalized Gibbs ensemble is the relevant ensemble for quenches with quadratic fermionic 
models, but a similar situation occurs for the usual statistical ensembles, i.e., microcanonical, 
canonical and grand-canonical. 

§ For the GGE ensemble things might be worked out by an appropriate representation of the Dirac’s 
delta, or by computing the moment generating function of Pgge(^.)- These tricks however would not 
work, for instance, for the microcanonical distribution. 
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proposed in 2001 by F. Wang and D.P. Landau, is a Monte Carlo method designed 
to compute the density of states of a classical statistical mechanics problem. The 
algorithm performs a non-Markovian random walk to build the density of states 
by overcoming the prohibitively long time scales typically encountered near phase 
transitions or at low temperatures. Besides the classical Ising and Potts models 
studied in the original papers usiEniizu, the method has been applied to the solution 
of numerical integrals [22], folding of proteins |23| and many other problems. 

Here is the plan of the paper. In Sec. [2] we show how the WLA can be extended to 
compute weighted density of states. In Sec. [3] we show how the weighted-WLA can be 
used to compute distributions related to quantum quenches with quadratic fermionic 
models. Finally, Sec. [4] contains a summary and future perspectives. 


2. Weighted Wang-Landau algorithm 


Let us consider a system with a discrete configuration space, where configurations can 
be labeled with an index a. Given a physical observable A, we define its weighted 
(coarse-grained) density of states: 


,(A) = y, w a 6 A A a 


( 7 ) 


with w a a positive weight. Here &AA a is a Kronecker delta, or, if the possible values 
of A a are too dense to keep them all, a suitable histogram-window-function coarse- 
graining of the Dirac’s delta. When w a = 1, we recover the usual density of states 
p(A), and the WLA can be used to estimate it 19] .2(1 [2T]. We will now show that, 
by properly modifying the WLA, we can compute p w (A) for generic w a s. 

To understand the gist of the approach, consider a generic positive function p w (A) 
which is our best guess for the desired p w (A) —, and set up a Markov chain random 
walk in which, given a state a, a new state a' is generated with a trial probability 
T(a'\a), which we will take to be symmetric, T(a' |a) = T(a|cd), and accepted with 
probability: 


R(a'\a) = Min 


1 Wa’ Pw{A a ) T(a\a') 
’ w a p w (A Q /) T(a'\a) 


( 8 ) 


With this standard Metropolis Monte Carlo prescription, we know that, after an 
initial transient, we will visit the configurations a with an equilibrium distribution 
P® q fulfilling the detailed balance condition and given by: 

pe q _ pi Wa 
a Pw(A«) ’ 

where C is a normalization constant. As in the WLA OS, while the random walk 
goes on, we collect a histogram h(A), updating h(A a ) —»• h{A a ) + 1 at each visited 
state a. At equilibrium, after N s steps, the “mean” histogram will then be given by: 


HA) = N s T.^ q SAA. 

a 




W 0 


Pw ( A a ) 


8 aa q 


= N S C 


Pw (A) 
Pw(A) ' 


( 9 ) 


Exactly as for the WLA ns, if our guess for p w (A) is a good approximation to p w (A), 
the histogram h{A) will be “almost flat” (see below). Obviously, during the random 
walk, together with the histogram h(A) we also update our guessed p w (A). Therefore, 
closely inspired by the WLA [19], we propose the following algorithm: 
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(0) Fix a modification factor / > 1, and set ln/5 w (A) = 0 and h(A) = 0 for all values 
of A; 

(1) Start the Monte Carlo procedure using Eq. <[8j) and update at each step the 
histogram and the weighted density of states with the rules h(A a ) —>• h{A a ) + 1 
and ln/5 w (A a ) ->• lnp w (A a ) + In/; 

(2) Stop the random walk when h(A) is “almost flat” (for instance [15] . when 
h(A) > 0.8 h for all values of A, where h is the mean histogram value). For the 
previous observations, at the end of this step In p w (A) is a good approximation to 
lnp w (A) with a discrepancy of order In/; 

(3) Reduce the value of / —> y/f, reset h(A) = 0 and restart the procedure from step 
(1) using the p w (A) just obtained. Stop this loop when In/ is smaller than the 
desired discrepancy e. 

A similar extension of the WLA has been already been introduced for the particular 
case in which w a is the Boltzmann distribution, with the aim of computing the free 
energy profile as a function of a reaction coordinate J23, (25] • In the present paper, we 
will use this algorithm to compute distributions where the weights are not Boltzmann- 
like, but rather associated to quantum quenches. 

Let us return for a moment to the original WLA. A first trivial observation is 
that, as it should be, the weighted-WLA with w a = 1/N coincides with the WLA. In 
many situations, when the size of the configuration space is too big and the density of 
states ranges over too many orders of magnitude, it is convenient, in computing p(A), 
to run many WLA over small domains = [A^j n , Anix ] • But then the update 
rule of the standard WLA has to be changed to avoid that, during the random walk, 
A a leaves the domain A^. This trick was already used in the first papers by Wang 
and Landau, when dealing with the largest sizes [15]. To avoid leaks from A^, the 
empirical solution was to reject any proposal to states a' with A a > £ A^, without 
any update of p{A) and h(A). With this prescription, however, there are “boundary 
effects”, actually a systematic underestimation of the density of states at the borders 
of the intervals [25]. Schulz et al. [25] showed phenomenologically that such boundary 
effects are eliminated by using the rather obvious update rule: given a proposal a', if 
A a i is outside the interval we remain in a and we update h.(A) and In p(A) using the 
state a, otherwise we accept a! with the usual rule. This update rule is just what is 
obtained, rigorously, by using our weighted-WLA. Indeed, the density of states in a 
restricted range A^ is proportional to a weighted density of states in which w a = 1 
when A a £ A^\ and zero otherwise. With these weights, the update rule of our 
weighted-WLA is exactly the one obtained phenomenologically by Schulz et al. [26]. 

3. Quantum quenches 

In this section we come back to the initial problem of computing the distributions 
p D {A) and p GGE {A) related to quantum quenches. We will show that with the 
weighted-WLA we can compute these distributions for sizes inaccessible with an 
exhaustive enumeration. 

We concentrate on quantum quenches in two models possessing a free fermionic 
description. The first model we considered is the fermionic Anderson model with 
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disorder in the local potential: 

L L 

H A = -t^2 (ctcj+i + h.c.) + h jc]cj , ( 10 ) 

t=i 1=1 

where cj (dj) creates (destroys) a fermion at site j, t is the nearest-neighbor hopping 
integral and hj is an uncorrelated on-site random potential uniformly distributed 
in the range [—W/2, W/2]. We assume periodic boundary conditions. It has been 
mathematically proven [27] that, for Hamiltonians like H A and in presence of any 
W > 0, all the single-particle eigenstates of H A are exponentially localized. The 
second model we considered still describes spinless fermions hopping on a chain, but 
now the hopping is long-ranged [25]: 

J^lrh = ^ ' tj 1 j 2 Cj 2 + h.c.) , (11) 

I1I2 

where tj 1 j 2 is a (real) hopping integral between sites j 1 and j 2 . We will take the tj 1 j 2 ’ s 
to be random and long-ranged, with a Gaussian distribution of zero mean, ( tj 1 j 2 ) = 0, 
and variance given by: 



Here 7 is a real positive parameter setting how fast the hoppings’ variance decays 
with distance. Notice that, for j\ = j 2 , we have = 1 for any 7 , hence the 

model has also on-site Gaussian disorder; by increasing the distance between the two 
sites \ji — j’ 2 1, the variance of the hopping integral decreases with a power law. The 
peculiarity of this long-range-hopping model is that, although being one-dimensional 
and regardless of the value of j3 (which hereafter is fixed to 1), it has an Anderson 
transition from (metallic) extended eigenstates, for 7 < 1 , to (insulating) power-law 
localized eigenstates for 7 > 1 E3E21EB. Physically, this is due to the fact that, 
for small 7 , long-range hoppings are capable of overcoming the localization due to 
disorder. Having access, in the same model, to physical situations in which the final 
eigenstates are extended (7 < 1 ) or localized (7 > 1 ) will clearly show the role that 
spatial localization plays in disrupting the ability of the GGE to describe the after- 
quench dynamics. Physically, spatial localization prevents the different “modes” of 
the system from having an infinite reservoir. 

For the considered quenches, we use as initial Hamiltonian Hq the clean chain 
with W = 0 and the same boundary conditions of the final Hamiltonian, i.e., periodic 
boundary conditions when quenching to H A and open boundary conditions when 
quenching to Hhp. The corresponding initial state |T 0 ) will be the filled Fermi sea, 
i.e., the ground state of H 0 with Np = L/ 2, where Np is the number of fermions. The 
reason behind this simple choice for Hq is that the “stationary state” reached does not 
depend, qualitatively, on the initial Hamiltonian being ordered or not, see Ref. m- 
The final Hamiltonian will be the Anderson model H A with W = 2, or the long-range 
hopping chain Hi r p with 7 = 0.5 or 2. In all cases the particle number is a constant of 
motion, therefore Np = L/ 2 for any time t > 0. To get a smoother size dependence 
of the computed quantities, the smaller size realizations are obtained by cutting an 
equal amount of sites at the two edges of the largest realization. 
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The two Hamiltonians, being quadratic in the fermion operators, can be 
diagonalized for any chain of size L in terms of new fermionic operators 

L 


4 = u oA ’ 

3=1 


(13) 


where Uj M are the wave functions of the eigenmodes of energy e^: H\/\rh = 

The energies and the associated wave functions Uj^ are obtained, for any given 
disorder realization of a chain of size L by numerically diagonalizing the L x L one- 
body hopping matrix. 

Given an observable A , consider the two distributions introduced before: 


Pd (A) =J2\c a \ 2 S{A-A aa ) 

a 


e ^ A ^ n 2 

Pgge(^) = / v y 5(A. A aa ) , 


(14) 

(15) 


where S(x) is the Dirac’s delta, {|a)} are the many-body eigenstates of H, A a p = 
{a\A\j3), c a = (al^o), and n“ = 0,1 is the occupation of the single-particle eigenstate 
p in the many-body eigenstate |a). These functions give the weighted distributions of 
A aa in the diagonal and GGE ensembles. 

Let us discuss a few technical details on the implementation we made, before 
discussing the physics emerging from our calculations. Notice that the sum over a 
is effectively restricted to the canonical Hilbert space TLn with a fixed number of 
particles N = Np in the diagonal ensemble, since c a = (a|’I'o) = 0 if N a ^ Np. 
No such restriction is in principle present in the GGE case, where the sum over a 
runs over the grand-canonical Hilbert space. By definition, the distributions are such 
that ( A) d = fAp D (A)dA and (A) GGE = J dAp GGE (A)A, where the integration is 
over the domain of A aa . As customary in any numerical finite-size study, one really 
needs to consider a coarse-grained version of these distributions, obtained by splitting 
the domain of A into small intervals of amplitude A. Such a coarse-grained 
distribution has exactly the form of a weighted density of states, see Eq. ([7]), with 
w a = |Cq, | 2 /A in the diagonal case, and w a = e~ /(A Z GGE ) in the GGE case. 

The configuration space {|a)} (i.e., the canonical Hilbert space TLn for the diagonal 
distribution and the full Hilbert space for the GGE) over which the two weighted 
distributions are defined is discrete and grows exponentially with the system size. 
The weighted-WLA is therefore the appropriate tool for the numerical computation 
of /O d (A) and p GGE (A). The eigenstates |a) which appear in the definition of p D (A) 
have a fixed number of fermions Np (the same of the initial state), while in / 0 GGE (A) 
the number of particles can change. In the weighted-WLA, for the diagonal ensemble, 
we use therefore a “particle conserving” proposal scheme: given a state |a), the state 
\a') is given by moving at random a fermion in one of the unoccupied single-particle 
eigenstates. In this case, the ratio w a < /w a which appears in Eq. ([ 8 ]), is equal to: 

_ |Cq '| 2 

| C q | 2 

where the coefficient |c a | 2 is the square of the determinant of a Np x Np matrix (see 
El App. D] for the explicit expression of |c a | 2 ). For the GGE case, instead, we do 
not have restrictions on the number of fermions and, given a state |a), we generate a 
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state |a') by changing the occupation of a randomly selected single-particle eigenstate 
p,. In this case: 

^ =e ±A, j 

W a 


where the + (—) sign appears when the mode p is initially occupied (empty). Let us 
recall that the Lagrange’s multipliers are obtained by requiring (\E' 0 |%^ai|'I'o) = 
(u^d/j) gge- This condition, written explicitly, reads: 


= 


- 1 


K t «] l 


( 16 ) 


where u° and u are L x L matrices whose elements u® v and Uj^ are the single-particle 
wavefunctions of the initial Hamiltonian Hq and the final one H, and = 0,1 is 
the occupation of the nth eigenstate of Hq in the initial state. The difference in the 
computational effort on computing the ratio w a >/w a in the two ensembles is evident: 
in the diagonal case at each step we have to compute the determinant of a Np x Np 
matrix, while in the GGE we have just to recover the value of e A ^ (they can be 
computed and stored before the Monte Carlo calculation because their number is L). 
Here we will show results for sizes up to L = 256, where both p D (/ 1) and p G GE (A) 
can be computed and compared. (For the GGE ensemble, we could reach L = 1024 
without problem.) In the numerical computations we used a minimum value of the 
WL parameter e = ln/ m j n = 10~ 6 , and we split the domain of A in L bins. Notice 
that the domain of A in p GGB (A) is always larger than the domain of p 0 {A) because, 
in the GGE, the many-body eigenstates do not have a restriction on the number Np 
of fermions. 

In the next two subsections we show the results obtained with the weighted- 
WLA for the calculation of p D (A) and p GGE {A) for two observables, the total energy 
and the local density. The physical picture emerging from the calculation of the full 
distribution function of the after-quench energy and local-density confirms and extends 
the results discussed in Ref. mm- In particular, we find clear differences between 
the diagonal and GGE distributions, even at the level of the variances, whenever a 
disorder-induced spatial localization is at play in the after-quench Hamiltonian. 


3.1. Probability distributions of the energy 

The first observable we consider is the total energy: Here A aa —» E a = 

where = (a|djj(i M |a) = 0,1 are the single-particle occupations of the eigenstate 
|a). In Fig. [l] we show the distributions In \p D (E)\/L and ln[p GGE (£’)]/L, computed 
for L = 128 and L = 256, for the three cases we have studied, i.e., quenches from an 
initially ordered half-filled chain Hq towards: 1) a long-range hopping Hamiltonian 
H\ r p with extended eigenstates (7 = 0.5, top), 2) H\ rh with localized eigenstates 
(7 = 2, center), and 3) an Anderson Hamiltonian Ha with a disorder width W = 2 
(bottom). Observe, first, that the distributions p D {E ) and p GG e (E) shown in Fig. [I] 
have identical average (denoted by a solid vertical line) 

(H) D = jdEp D (E)E = JdEp G G E (E)E = (H) G Gv 

This result comes directly from the fact that the energy does not fluctuate in time 
(i.e., the diagonal energy coincides with the average energy (\E , 0 |-Hj’I'o)) and GGE 
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Figure 1. Value of \n[pr>(E)\/L and Iii[pgge(E)]/ L computed with the weighted- 
WLA. The gray curves are obtained with L = 128, while the black ones with 
L = 256. The solid vertical lines are the average energy after the quench, i.e. 

f° r L — 256. The three panels are obtained starting from the ground 
states of clean chains and quenching to different disordered Hamiltonians: panel 
(a) long-range hopping with 7 = 0.5 (extended eigenstates), panel (b) long-range 
hopping with 7 = 2 (localized eigenstates) and panel (c) Anderson model with 
W = 2. For the computations we used a single disorder realization and, to get a 
smoother size dependence, the smaller size realization is obtained by cutting an 
equal amount of sites at the two edges of the larger realization. These distributions 
are obtained for a single realization of the couplings, but we verified that, for large 
sizes, the results are self-averaging. 


fixes the occupation of the fermionic eigenstates in such a way as to exactly reproduce 
(^ol-ffl^o). The form of the two distributions, however, differs considerably, most 
notably at the extremes of the spectrum, and for the Anderson model case. Let us 
now consider the fluctuations of the energies in both distributions. In the diagonal 
ensemble the variance is: 

4,d = J dEp D (E) E 2 - (H)l = (H 2 ) d - (. H)l , (17) 

where the expression on the right-hand side holds only for the Hamiltonian (it would 
not apply to arbitrary operators, because ( A aa ) 2 ^ {a\A 2 \a)). An entirely similar 
expression applies to the GGE case. Since the energy is an extensive operator, it is 
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reasonable to ask what happens to the fluctuations in the energy-per-site e = E/L , 
which are simply given by a 2 eu = a 2 ED /L 2 , and (r 2 eGGE = o 2 Egg JL 2 . On pretty general 
grounds, for quenches of local non-integrable Hamiltonians, it is known DU IDl that 
a 2 0 —> 0 in the thermodynamic limit, L —> oo. Indeed, as shown in Fig. [ 2 ] both a 2 0 
and a 2 GGE decrease to 0 for L —)• 00 for the three considered cases. For our quadratic 
problems, however, we can say a bit more. First of all, from the explicit expression in 
Eq. @ after very simple algebra (mainly using Wick’s theorem), we arrive at: 

ff e,GGE = J2^Z e l n l 9 “ n °) ’ ( 18 ) 


a 


2 

e,D 


= a 


2 

e,GGE 


e A»i e M2 

#*17^*2 


I c 

| LT MlM2 


2 


(19) 


where G° lM2 = l^o) is the t = 0 one-body Green’s function. The off- 

diagonal elements of G ° lM2 play here an important role, and the second term in a 2 D 
originates from the fact that, by definition, GGE does not include correlations between 
different eigen-modes, i.e., = 0 , when /Hi 7 ^ ^ 2 - 


10 ° 


10 


,-l 


10 


,-2 


CN cu 

b 


10 


,-3 


10 


10 - 


,-4 


^ //(,-> H,rh (7 = 0.5, ext.) 

S'', 

''ft, 

''-a. 



■ , O. 
'*■. 


■ 


' ' 7>- , H 0 -4 H kh (7 = 2, loc.)' -ft.. 


Ho -> H a (W=2) 


0 „ 

“ A „ ' 


XX ^ 


" 10 - . 




10 1 


10 " 


10 3 


L 


Figure 2. Variances cr^ D = cr^ D / L 2 (empty circles) and cr^ GGE = < 7 ^ GGE /-^ 2 
(solid triangles) as a function of the size L. The data are obtained using the same 
set of quenches used in Fig. [I] and the values are computed using Eqs. ( |19[ ). Error 
bars are calculated by averaging over 20 different realizations of the disorder. 
The dashed lines are power law fits ~ L~ s , where s « 1 for the Anderson 
case, while, for the quench to H\ r h, s ~ 0.82 when 7 = 0.5, and s ~ 0.95 when 
7 = 2. Notice the observable difference between D and cr^ GGE when the final 
eigenstates are localized. 

Let us first consider the Anderson model case. Assuming, as done so far, a 
bounded distribution of disorder, we are guaranteed that a finite bound e max exists 
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such that |e M | < e n 


for any L. With this assumption, it is easy show that ct/ GGE has 


to go to zero at least as 1/L for L —> oo. Indeed, the occupation factors appearing in 
°e,GGE are suc ^ that 0 < n° (l - n°) < 1/4. Hence: 


a e,GGE — 


L 2 




(!-<)< 


4 L 


( 20 ) 


The same statement can be made for tTg D , because the difference between the two 
variances has a similar upper bound: 


-o£ggJ< ^2 E I 


G 


I — 


L 2 


E 


a 


MlM2 I — 


( 21 ) 


where we used that Y 

A2 


PlM2 


Ir 0 

M1M2 I 


T iv 


and 


= Np < L. Nevertheless, although both a 2 
e GGE go to 0 as 1/L for the Anderson model, they do so with a different pre-factor, 
see Fig. [2] and comments below. 

For the quenches to H\, -h, a bound e max for the single-particle spectrum is in 
principle not defined: one can think of rare realizations in which the hopping is large at 
arbitrarily large distances, which would give an unbounded distribution of eigenvalues 


e^. Indeed, the behavior of both cr^ 


and 


e,GGE 


if, 


that the power- 
(we find s ss 0.82 


suggest, see Fig. 

law approach to 0 might be slower than 1/L, i.e., as L~ s with s < 
for the case 7 = 0.5 and s ~ 0.95 for 7 = 2). While this might be a finite-size artifact, 
we find it intriguing that such deviations are quite clearly seen for quenches to H\ T h' 
they might be due to the power-law nature of the hopping integral variance. 


Concerning the similarity between crf D and er 


e,GGE> 


we observe that the two 
essentially coincide for the case of a quench to H\ v p with extended eigenstates, while 
there is a small discrepancy for the quench to H \ r h with localized eigenstates, and 


a quite clear different pre-factor in the Anderson model case, a/ 


C d /L and 


G n 


^e,GGE Cg G e/L with C GGE < 

analyzing the term e in e H 2 |G° lA(2 | 2 which appears in Eq. (19). In Fig. [ 3 J panel 

(b), we show the structure of the matrix |G° lAl , | 2 for the three quench cases. We divide 
this matrix into four sectors, one for each sign of the single-particle energies and 
e^ 2 : in two of these quadrants the product is positive (top-right and bottom- 

left), and in the others is negative. For quenches to H\ r p this matrix is almost equally 


This different pre-factor can be understood by 
2 


distributed in all the four sectors: the sum Y 
leading to <r 


| 2 has cancellations, 


matrix 

2 

°e,GGE 


IGJU 


e,GGE 

2 


M 1 AM 2 C /U C M2 I V ^ r /2l/22 

< 7 / D for large sizes. For quenches to Ha, on the contrary, the 


is mainly concentrated in the sectors in which 


< 0, leading to 


Finally, let us comment on one aspect of the distributions shown in Fig. [l] which 
can be easily understood from the single-particle occupations shown in Fig. [3} We see 
that, when the after-quench Hamiltonian is the Anderson model, Pn{E) has both mode 
(i.e., maximum value) and average very close to the ground state energy: the quench 
excites mostly the low-energy part of the many-body spectrum. On the contrary, for 
both the quenches towards F/i r h, mode and average are almost in the middle of the 
many-body spectrum; there, indeed, the quench is more dramatic: we are going from 
the ground state of a chain with nearest-neighbor hopping to a disordered chain with 
long-range hopping. This is evident by looking at the occupations = ('F 0 |dld^l'kg) 
as a function of the single-particle energy e^, shown in Fig. [ 3 } panel (a). By definition, 
only the eigenstates of H 0 with e/ < 0 are occupied in |*Fo). The quench to Ha only 
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Figure 3. Panel (a): occupations n? = |'I'o) as a function of the 

single-particle energy e: /y . Panel (b): representation of the matrix |G®„ 2 | 2 . For 
the diagonal and off-diagonal elements we add a black pixel when the value exceeds 
their mean value. For the diagonal elements the mean value is x = in® ) 2 /L, 
while for the off-diagonal elements the mean value is (ATp — xL)/L(L — 1), where 
jVp is the number of fermions in the initial state, and we used the relation 
£, 1M2 |g° 1M2 | 2 = (see |31l App. D]). The vertical and horizontal lines 
indicate the indexes at which the single-particle energies e Ml and e^ 2 change sign, 
and the signs shown in the four quadrants are those of the product e fJ , 1 e M2 . For 
the two panels we used L = 256 and the same quenches used in Fig. [I] and Fig-0 


slightly modifies the initial occupations: rN. apart for fluctuations due to disorder, goes 
smoothly from 1, in the lower part of the single-particle spectrum, to 0, in the highest 
part of the spectrum. On the contrary, for the quenches towards Hhh, the single¬ 
particle spectrum is entirely excited, both the positive and the negative energy part. 
This explains why, for these quenches, the after-quench energy ('F 0 |.fi|'I' 0 ) = XX e n n< ji 
is near the center of the many-body spectrum. 

3.2. Probability distributions of the local density 

Let us now consider the local density hj = CjCj, perhaps the simplest one-body 
observable. For definiteness, we concentrate on j = L/ 2, the center of the chain. 
It is important to stress that we are going to consider the fluctuations of hj before any 
possible average over the sites j: averaging over the sites j an intensive local operator 
would effectively send to zero the fluctuations in the thermodynamic limit | |13j , while 
we will show that, for a fixed j, finite fluctuations survive in the thermodynamic limit 
when the eigenstates are localized, due to disorder. 

The diagonal and GGE distributions p D {n) and p GGE (n) are now constructed using 
the matrix elements n aa = ( a\hj\a ) = \ u ju\ 2n u^ w ^ ere ^2 = 0,1 are J as before, 
the single-particle occupations of the eigenstate |a). In Fig. 0 we plot ln[p D (n)] and 
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ln[pGG K (n)] ; computed for the three quenches discussed before. The case of a quench 
to -ffuh with 7 = 2.0 (localized eigenstates) is quite peculiar. The values that n can 
assume is actually split in two separated domains, one just above n = 0 and one just 
below n = 1 , and the mean value is exactly in the middle, where no values of n aa 
happen to fall. This is due to the strong spatial localization of the eigenstates. As we 
show in Fig. |5] at fixed j, the value of \uj^\ 2 is strongly localized in a single eigenstate 
fb. This implies that the value n aa = V has a strong jump when we move 

from a state |a) in which = 0, to the state |a) in which n~ = 1. For the quench 
to Ha, with W = 2, the localization is not strong enough to produce such a gap: we 
however expect this to happen for larger values of the disorder amplitude W. 



Figure 4. Distributions of the local density hj at the center of the chain, 
j = L/2, for the diagonal ensemble and the GGE, and for L = 256 (black curves) 
or L = 128 (gray curves). In panel (a), we plot ln[pD(n)]/L and ln[pGGE(^)]/T, 
while in panels (b) and (c) ln[pD(n)] and ln[pGGE (?"*•)]• The vertical lines are the 
diagonal and GGE average of hj, which coincide for the local density. The three 
panels are obtained using the same quenches of Fig. [T] 

Since hj is a one-body operator, the diagonal and GGE averages coincide El, 
and therefore, the mean value of the two distributions is the same: 


dn p D (n) n = 


dnp GGE (n)n . 
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Figure 5. Squared single-particle wavefunction |itj M | 2 as a function of the 
eigenstates index /.t, at fixed site j = L/2. We have taken L = 256 and the 
three panels are obtained using the same quenches of Fig. [T] 

We also note that, for this observable, hj 1 = fij for any positive integer to, and 
therefore (n™) D = (^™)gge- However, this does not allow us to conclude that the two 
distributions p D (n) and p GGE (n) coincide, since, unlike the case of the total energy, we 
have that, for instance: 

Jdnp D / gge (iT) n 7^ (hy) d/gge • 

The variance of the two distributions can be computed by exploiting again Wick’s 
theorem. We find that: 

^n.GGE = K> I 4 n ° ( X - n V) ( 22 ) 

A 4 

C r n,D=^n,GGE- l U J>i| 2 l U J> 2 | 2 |G> lM2 f' - ( 23 ) 

In Fig. [ 6 ] we plot cr^ GGB and cr^ D as a function of size. We see that, in both 
ensembles, the variances vanish as 1/L when quenching to i?i r h with 7 = 0.5 (extended 
eigenstates) while they are finite when quenching to Ha and to H] r h with 7 = 2 , i.e., 
when the final Hamiltonian has localized eigenstates. These results agree with the 
findings of Ref. [52] . who show that, for large L, the variance of few-body intensive 
(but not site-averaged) observables remains finite both in the microcanonical ensemble 
and in the diagonal ensemble for the Aubry-Andre model. 
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From the equation for cr^ GGE , we see that it is related to an inverse participation 
ratio (IPR): the sum is over the eigenstates /z, each /z weighted with the corresponding 
occupation factor 0 < (l — n°) <1/4 depending on the initial state. It is therefore 
easy to realize that: 


f 7 n,GGE — ^ ^ ' 


IPR, 




(24) 


where the last equality defines the IPR at fixed site j. This shows that, whenever the 
IPR goes to zero, i.e., when the final Hamiltonian has delocalized eigenstates, cr^ GGE 
goes to zero as well. For a final Hamiltonian with localized eigenstates we have instead 
the opposite: there is at least one eigenstate /z localized around j, and therefore there 
is a single-particle wavefunction Ujp, which does not vanish in the thermodynamic 
limit; if the initial occupation rC- of this eigenstate is such that 0 < n? < 1, then 
a n gge remain finite in the thermodynamic limit. 

Concerning <r^ D , Eq. (231 can be rewritten as: 


' n,D i 
2 


_ 2 _ s:2 

a n, GGE °jj 


(25) 

where <5^ denotes the mean squared time-fluctuations of the single-particle Green’s 
function G JU - 2 (t) [ID] : 


nj2 


= lim - 

t—too t 


d a \SG jlh (t')\ 


(26) 


dGj 1 j 2 (t) = Gj 1 j 2 (t) — Gj 1 j 2 being the time fluctuation with respect to the long-time 
average Gj, j2 . Physically, <$? ■ is the averaged long-time fluctuation of the local density 

hj = Cjdj. In Refs. nni nu we have shown that if the final Hamiltonian has extended 
eigenstates, then dC « 1 /L for large sizes, while 6jj remains finite when the final 
Hamiltonian has localized eigenstates. This explains all the features shown in Fig. [6} 
in particular the clear difference between D and < 7 ^ GGE in all cases. 


4. Summary and conclusions 

In this paper we have introduced a Monte Carlo method obtained by a rather 
natural extension of the Wang-Landau algorithm US HB HU which allows to 
compute quite general weighted distribution functions of the form relevant to quantum 
quenches, see Eq. |6]). We have used this approach to analyze quantum quenches 
for free-fermion Hamiltonians in presence of disorder. For these systems, thanks to 
Wick’s theorem, after-quench expectation values and time averages require a modest 
computational effort, proportional to a power-law of the size L HU. However, the 
calculation of full probability distributions — like the diagonal ensemble distribution 
p D (H), Eq. @, or the GGE one pgge(41), Eq. ([5]) — would still require a sum over 
an exponential number of terms, hence unfeasible beyond very small sizes. 

Although quadratic, hence with an extensive number of conserved quantities, 
these free-fermion problems are not described by the GGE ensemble whenever the 
disorder is such that the after-quench eigenstates are localized. More precisely, while 
the GGE ensemble is known to correctly capture the long-time average of any one- 
body operator, almost “by construction” HU , it does not capture correlations induced 
by the spatial localization of the eigenstates. Our study further explored this issue 
by explicitly calculating and comparing the full probability distributions of both the 
energy and the local density in the two relevant ensembles. 
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Figure 6. Plot of the variances <x^ D (circles) and n 2 n GGE (triangles) as a function 
of size. The data are obtained using the same set of quenches used in Fig. [2] and 
the values are computed using Eq. l ]23[ i. Error bars are calculated by averaging 
over 20 different realizations of the disorder. The dashed lines are power-law fits 
rrj) ~ L~ s , where s ~ I iri both cases. 


Concerning the energy, we have explicitly verified that the form of the two 
distributions for the diagonal and GGE ensembles differs considerably, most notably 
at the extremes of the spectrum, and for the Anderson model case. More in detail, 
we have verified that, regardless of the final Hamiltonian, the averaged fluctuations of 
the energy-per-site, [cr^ D ] av and [ffe GGE ] av ' go to zero in the thermodynamic limit, 
see Fig. [2j in agreement with the general analysis of Refs. mm- Nevertheless, we 
find that there is a clearly detectable difference in the two variances when the final 
Hamiltonian has localized eigenstates. 

In addition to the energy, we studied the local density distributions. For this 
observable, it was already known that, even in presence of disorder and localization, 
the GGE expectation value coincides with the diagonal average m, a property true, 
more generally, for any one-body operator m- Our numerical results confirm that 
even if the averages of pu(n) and pgge(«) coincide, the two distributions are different 
when localization is present, with clearly detectable differences already at the level 
of the variance, see Fig. |6j cr^ GGB and cr^ D differ by a quantity which represents 
the averaged long-time fluctuations of the local density maun, which remain finite 
whenever the final Hamiltonian has localized eigenstates. 

Other many-body operators, like density-density correlations, might be analyzed 
in a similar way. Here, even the average values are not in general well described by the 
GGE distribution whenever localization is at play m we expect, once again, clearly 
visible discrepancies between the diagonal and GGE distributions in such cases. 

In conclusion, as explained above, the weighted-WLA we have presented 
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circumvents the difficulty associated to the exponentially large Hilbert space to be 
visited, even when the relevant ingredients entering the distribution — matrix elements 
and overlap between states — can be calculated quite effectively. In principle, the 
applicability of the method is not limited to “quadratic fermion problems” of the 
type we have considered in the paper: If, by Bethe Ansatz, or any other exact 
technique or even by a suitable quantum Monte Carlo approach, one would be able to 
calculate matrix elements and overlaps, the method we have illustrated would provide 
an effective Monte Carlo sampling of the relevant distribution functions. 
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